
###############################################################################
### Did 3G Make Afghan Insurgents Fight More Effectively? 
### Replication files for analyses
### Mehmet Erdem Arslan
###
### Appendix E: Heterogeneous Effects  
###############################################################################

# setwd("your file path here")

data <- readRDS("R/Heterogeneous_effects.rds")

# id: Admin2 (district id)
# time: quarter
# all: number of insurgent attacks (source: Global Terrorism Database - GTD)
# log_all: log number of insurgent attacks (source: GTD)  
# crd: number of coordinated attacks (source: GTD)
# two: proportion of the area covered by 2G network
# three: proportion of the area covered by 3G network
# lag.two: proportion of the area covered by 2G network, lagged for one quarter
# lag.three: proportion of the area covered by 3G network, lagged for one quarter
# POP_2000: population 
# log_pop: (log) population 
# DIST2PROVCAP: distance to province capital 
# nl_mean: nighttime lights 
# ROAD_DENSITY: road density 
# lag_all: number of insurgent attacks, lagged for one quarter
# laglag_all: number of insurgent attacks, lagged for two quarters 
# *For information on population, night lights, distance to province capital, and road density, data from Zhukov et al. (2019) are used. Population data is for the year 2000, and population, distance to province capital, and road density variables are time-invariant.

library(fixest)
library(marginaleffects)
library(ggplot2)
library(texreg)
library(grid)
library(gridExtra)

# (log) population 
twfe_pop <- feols(log_all ~ lag.two*log_pop + lag.three*log_pop | id + time, data = data)

# distance to province capital
twfe_dist <- feols(log_all ~ lag.two*DIST2PROVCAP + lag.three*DIST2PROVCAP | id + time, data = data)

# night-lights 
twfe_nl <- feols(log_all ~ lag.two*nl_mean + lag.three*nl_mean | id + time, data = data)

# road density 
twfe_road <- feols(log_all ~ lag.two*ROAD_DENSITY + lag.three*ROAD_DENSITY | id + time, data = data)

# past insurgent violence (t-1)
time_lag_dv <- feols(log_all ~ lag.two*lag_all + lag.three*lag_all | id , data = data)

# past insurgent violence (t-2)
time_laglag_dv <- feols(log_all ~ lag.two*laglag_all + lag.three*laglag_all | id , data = data)


#### Plots for heterogeneous effects ####

plot1 <- plot_predictions(twfe_pop, condition = list("lag.three", 
                                                     log_pop = c(round(quantile(data$log_pop), 2))), 
                          conf_level = 0.9, rug = TRUE) + ylim(0, NA) + 
  theme_classic() + 
  scale_x_continuous("3G coverage") + 
  scale_y_continuous("(log) Number of attacks") +  
  guides(colour = guide_legend(title = "(log) Population"), 
         fill = guide_legend(title = "(log) Population"))


plot2 <- plot_predictions(twfe_nl, condition = list("lag.three", 
                                                    nl_mean = c(round(quantile(data$nl_mean), 2))), 
                          conf_level = 0.9, rug = TRUE) + ylim(0, NA) + 
  theme_classic() + 
  scale_x_continuous("3G coverage") + 
  scale_y_continuous("(log) Number of attacks") +  
  guides(colour = guide_legend(title = "Night lights\n(mean)"), 
         fill = guide_legend(title = "Night lights\n(mean)"))

plot3 <- plot_predictions(twfe_dist, condition = list("lag.three", 
                                                      DIST2PROVCAP = c(0.61, 28.55, 52.35, 79.64, 200)), 
                          # > round(quantile(data$DIST2PROVCAP), 2)
                          # 0%    25%    50%    75%   100% 
                          # 0.61  28.55  52.35  79.64 307.29 (max is one outlier, Wakhan, Badakhshan)
                          conf_level = 0.9, rug = TRUE) + ylim(0, NA) + 
  theme_classic() + 
  scale_x_continuous("3G coverage") + 
  scale_y_continuous("(log) Number of attacks") +  
  guides(colour = guide_legend(title = "Dist. to Province\nCapital"), 
         fill = guide_legend(title = "Dist. to Province\nCapital"))

plot4 <- plot_predictions(twfe_road, condition = list("lag.three", 
                                                      ROAD_DENSITY = c(round(quantile(data$ROAD_DENSITY), 2))), 
                          conf_level = 0.9, rug = TRUE) + ylim(0, NA) + 
  theme_classic() + 
  scale_x_continuous("3G coverage") + 
  scale_y_continuous("(log) Number of attacks") +  
  guides(colour = guide_legend(title = "Road density"), 
         fill = guide_legend(title = "Road density"))


plot5 <- plot_predictions(time_lag_dv, condition = list("lag.three", 
                                                        lag_all = c(0, 5, 10, 15 ,20)), 
                          conf_level = 0.9, rug = TRUE) + ylim(0, NA) + 
  theme_classic() + 
  scale_x_continuous("3G coverage") + 
  scale_y_continuous("(log) Number of attacks") +  
  guides(colour = guide_legend(title = "N. of attacks\n at t-1"), 
         fill = guide_legend(title = "N. of attacks\n at t-1"))

plot6 <- plot_predictions(time_laglag_dv, condition = list("lag.three", 
                                                           laglag_all = c(0, 5, 10, 15 ,20)), 
                          conf_level = 0.9, rug = TRUE) + ylim(0, NA) + 
  theme_classic() + 
  scale_x_continuous("3G coverage") + 
  scale_y_continuous("(log) Number of attacks") +  
  guides(colour = guide_legend(title = "N. of attacks\n at t-2"), 
         fill = guide_legend(title = "N. of attacks\n at t-2"))


#### Figure 6 ####
pdf("Outputs/Figure6.pdf", height = 8, width = 9.5)
grid.arrange(plot1, plot2, plot3, plot4, nrow = 2, ncol = 2)
dev.off()

#### Figure 7 ####
pdf("Outputs/Figure7.pdf", height = 4, width = 9.5)
grid.arrange(plot5, plot6, nrow = 1, ncol = 2)
dev.off()

#### Table 11 ####
texreg(list(twfe_pop, twfe_nl, twfe_dist, twfe_road, time_lag_dv, time_laglag_dv), 
       single.row = FALSE,custom.header = list("DV: All Insurgent Attacks (GTD)" = 1:6), 
       digits = 3, booktabs = TRUE, dcolumn = TRUE, scalebox = 0.8,  
       custom.model.names = c("Population", "Night Lights", "Dist. to Prov.Cap.", 
                              "Road Density", "Ins. Viol. (t-1)", "Ins. Viol. (t-2)"))


rm(list = ls())

### end of script ###
